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Fluctuations and noise may alter the behavior of dynamical systems considerably. For example, 
oscillations may be sustained by demographic fluctuations in biological systems where a stable fixed 
point is found in the absence of noise. We here extend the theoretical analysis of such stochastic 
effects to models which have a limit cycle for some range of the model parameters. We formulate a 
description of fluctuations about the periodic orbit which allows the relation between the stochastic 
oscillations in the fixed point phase and the oscillations in the limit cycle phase to be elucidated. 
In the case of the limit cycle, a suitable transformation into a co-moving frame allow fluctuations 
transverse and longitudinal with respect to the limit cycle to be effectively decoupled. While lon- 
gitudinal fluctuations are of a diffusive nature, those in the transverse direction follow a stochastic 
path more akin to an Ornstein-Uhlenbeck process. Their power spectrum is computed analytically 
within a van Kampen expansion in the inverse system size. This is carried out in two different ways, 
and the subsequent comparison with numerical simulations illustrates the effects that can occur due 
to diffusion in the longitudinal direction. 

PACS numbers: 05.40.-a, 02.50.Ey, 82.40.Bj 



I. INTRODUCTION 

The effect of noise on nonlinear dynamical sys- 
tems has been studied for some time jl| and is now a 
substantial field, yet significant new aspects continue 
to be unearthed. One of the most recent of these 
concerns systems which fundamentally involve dis- 
crete entities, for example individuals in an ecological 
system. These have populations which are modeled 
stochastically, for example random births and deaths. 
In such cases it may be that stochastic effects alter 
the behavior of non-linear systems substantially, and 
crucial differences between the properties of a given 
system can be observed in the presence and in the 
absence of noise. Examples can be found in the con- 
text of predator- prey population dynamics 0, @| , in 
evolutionary game theory [1, 0], in cyclic trapping 
reactions [oL in models of opinion dynamics Q, in 
epidemics [g, [^, [lO'] and in connection with genetic 
networks or biochemical clocks [H, [l^ . In these 
cases the "noise" is intrinsic to the system itself; in 
the parlance of ecology it is 'demographic stochastic- 
ity', rather than environmental stochasticity. 

One of the most intriguing effects found in these 
systems concerns the existence of oscillatory behav- 
ior. It has long been conjectured that in some sit- 
uations the influence of noise due to demographic 
stochasticity would be sufficient to perturb the sta- 
tionary state, predicted by a deterministic or mean- 
field type analysis, to produce cyclic behavior (T3 |. 
Oscillatory behaviors of this kind are referred to as 
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quasi-cycles [15j. This effect can be demonstrated in 
a simple and straightforward way, and analytic re- 
sults derived show very good agreement with sim- 
ulations This latter study conclusively demon- 
strates that while the deterministic limits of models 
of these systems exhibit stable fixed-point behavior, 
their stochastic analogues show resonant and persis- 
tent oscillations about these fixed point solutions. 

The starting point for formulating models of this 
kind is to define the system with a finite number of 
constituents (e.g. molecules in a chemical reaction 
system, individuals in the context of population dy- 
namics, or agents in models of social dynamics) which 
interact according to a given set of possible reactions, 
whose occurrence is determined by random factors 
(see e.g. [l^ [l^ and references therein). In spa- 
tial diffusive systems, for example, a certain molec- 
ular reaction may occur only if all necessary reac- 
tants are present at a given site in space and time, 
and similarly a predator in a model of population dy- 
namics may feed upon a unit of prey, only if both 
meet. Events such as birth and death typically oc- 
cur at random with Poisson statistics in such sys - 
tems, providing another source of stochasticity [17| . 
This is demographic stochasticity. In the limit of in- 
finite particle numbers, such systems are faithfully 
described by deterministic equations, also called rate 
equations. These ordinary differential equations for 
concentrations of the different reactants address the 
behavior of the system on a mean-field level. As- 
suming a well-mixed population the rate equations 
are zero-dimensional and describe uniform densities 
as functions of time. The stochasticity present on 
the level of interactions between individuals is aver- 
aged out in this case of an infinite system size. On 
the other hand, a systematic study of first-order cor- 
rections to the rate equations due to finite system 
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size can be made. Such studies have captured the 
behavior of fluctuations about the stationary mean- 
field solution which are small enough that a linear 
approximation is sufficient. In several examples, the 
power spectra of the fluctuations have been computed 
analytically [2,. JS]. 

The aim of the present paper is to extend these 
existing analytical descriptions of flnite-size stochas- 
tic effects to systems which on the mean-field level 
do not always approach a stable fixed point, but in- 
stead may tend towards a stable and periodic limit 
cycle solution for some range of the parameters. The 
effects of demographic stochasticity on such systems 
have been studied numerically, for example in ^ for 
predator-prey systems with a non-linear functional 
response. There the shape of the resulting auto- 
correlation functions of predator or prey densities 
have been used to help to distinguish between noisy 
limit cycles and quasi-cycles. This problem could 
equivalently be analyzed by looking at the power 
spectra. Distinguishing these two types of cyclic mo- 
tion might also be of relevance in the context of bio- 
chemical clocks and genetic networks 11]. Also of 
interest is to understand what happens at the bound- 
ary between the regimes where there is a stable fixed 
point and where there is a stable limit cycle. Do the 
cycles continuously merge with each other or are they 
unrelated? 

To develop the analytical tools required to study 
such systems, we focus on the Brusselator model 
[H [13, [H, 111, HH, Hi]. This model is a sim- 
ple example of an auto-catalytic, oscillating chemical 
reaction [isl . [23j . Auto-catalytic reactions are those 
in which the presence of a given reactant acts to in- 
crease the rate of its own production. A real-world 
realization of oscillatory chemical reactions is given 
by the celebrated Belousov-Zhabotinsky reaction [IJ] • 
The corresponding deterministic rate equations are 
known to have limit cycle solutions provided that 
model par ameters (i.e. reaction rates) are suitably 
chosen (l6l [2ll | . It also has the properties required to 
show resonant oscillations in the regime where a sta- 
ble fixed point exists. Therefore it has the necessary 
features required for our investigation. 

Earlier studies of the Brusselator system, e.g. by 
Tomita et al. |22| or by Scott et al. [ll| have ad- 
dressed finite-size corrections to the dynamics of this 
system based on van Kampen expansions of the corre- 
sponding master equation [25,] , but to the best of our 
knowledge no systematic attempts have been made 
to study temporal correlations (i.e. autocorrelation 
functions or power spectra) of this system. Specif- 
ically, our objective is to build on the calculations 
performed in [22|, and in particular to address these 
temporal correlations analytically. We also present 
a systematic account of different schemes of trans- 
formation into a co-moving frame (following the mo- 
tion of the deterministic system around the limit cy- 



cle). We link these with different interpretations of 
stochastic simulations based on the Gillespie algo- 
rithm [1^. It should be noted that recent studies 
by Frey et al. 0, [| also address models with limit 
cycles but those expansions are performed about un- 
stable fixed points in the interior of the limit cycle, 
whereas an expansion about the cycle itself is carried 
out in the present work. 

The remainder of the paper is organized along the 
following lines. In section |TT1 we introduce the model 
system initially as a mean field model and later as 
an individual based model. We then introduce the 
techniques of the van Kampen system-size expansion 
in Section Ulll bv applying them to the case where the 
mean-field dynamics approach a fixed point. A lin- 
earization about a stable limit cycle solution is carried 
out in Section IIVI we describe the limit cycle itself, 
the analysis of the Floquet multipliers and exponents 
and finally we discuss the co-moving Frenet frame. In 
Section|V]we use these tools to study the full stochas- 
tic problem of large but finite systems containing a 
limit cycle in the mean field. We derive the power 
spectrum describing the fluctuations in this case and 
also compare this prediction against numerical simu- 
lations of the individual-based Brusselator model un- 
der two different possible interpretations. In the final 
section we summarize our results and provide an out- 
look for future work. 



II. THE BRUSSELATOR MODEL 

A. Deterministic mean-field rate equations 

The non-spatial Brusselator model is well known 
in the form of the mass-action kinetic equations in 
two dynamic variables xi{t) > and X2(t) > repre- 
senting the time-dependent concentrations of a pair 
of d yna mic reagents. These equations are of the form 

ii = 1 - xi{l + b ~ CX1X2), 

±2 = Xi{b - CX1X2), ^ ' 

where h and c are constant and positive model pa- 
rameters related to the reaction rates. Since they will 
form the basis of the work we will describe in this pa- 
per, we will briefly review the structure of these equa- 
tions. Further details may be found in textbooks on 
nonlinear dynamics, for instance [2l| . For simplicity 
we will use the shorthand x(t) = {xi{t) , X2{t)) for the 
two-dimensional vector of concentrations. 

The two-dimensional and first-order nature of the 
Brusselator restricts the class of solutions to three 
possibilities. We may have fixed points at which 
the left-hand side of Eq. ((T|) vanishes, limit cycles 
where the dynamics repeats in a periodic orbit or 
a unbounded behavior where the solution tends to- 
wards infinity and never returns. It turns out that 
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this last possibility can be eliminated for the Brus- 
selator if we require that the dynamics begins within 
the positive quadrant (xi > and X2 > 0) and also 
that we have positive and finite parameter values. A 
limit cycle does exist for a range of parameter values; 
this will be discussed in Section IIV Al This leaves 
fixed points, which we will discuss here. It turns out 
that there is only one fixed point of Eq. ([T]) given by 
X* = {xl,X2) = (l,(6/c)). The behavior of the sys- 
tem close to this fixed point can be understood within 
a linear stability analysis. To this end, one considers 
a deviation e^{t) = x(i) — x* from the fixed point. 
The pre-factor e here indicates that we assume these 
deviations to be small. Later, in the context of a van 
Kampen expansion, this parameter will take on a spe- 
cific meaning in terms of the system size. Assuming 
that an expansion to linear order is appropriate, one 
then has 

i = K*i, (2) 
where the Jacobian at the fixed point is given by 

The eigenvalues of the matrix K* determine the sta- 
bility or otherwise of the fixed point {xl,X2)- They 
are both found to be real so long as |6 — 1 — c| > 2^/c. 
Otherwise the eigenvalues form a complex conjugate 
pair with a real part which is negative if 6 < 1 + c. 
The line b = 1 -I- c is a family of Hopf bifurcations 
which separate the parameter space into two phases: 
one in which there exists a single globally stable fixed 
point (b < 1+c) and another in which there is a single 
globally stable limit cycle (6 > 1 + c). The resulting 
phase diagram is depicted in Fig. [T] 

B. Microscopic multi-particle dynamics 

We can also discuss the Brusselator on the level of 
individual molecules. In this case, the system is de- 
scribed by four chemical reactions between four sub- 
stances, A,B,Xi and X2 [H, [U, HI . The number 
of A and B particles is, by construction, constant in 
time, so that their populations do not form degrees 
of freedom; their role is merely to set the reaction 
rates. The state of the system at any time is there- 
fore described by the number of molecules of each 
of the two substances Xi and X2, denoted by ni{t) 
and 77,2 (i). Both of these are non- negative integers 
at any time and we will define n = (ni,n2). The 
molecules are considered to interact randomly due to 
short-timescale fiuctuations such as thermal activa- 
tion. Hence we model their occurrence with Pois- 
son statistics. As such we need only specify the ex- 
pected number of occurrences per unit time for each 
reaction. We denote these transition rates as for 
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FIG. 1: (Color on-line) Phase diagram of the Brusselator 
fixed point at x* — (1, ib/c)) as obtained from the de- 
terministic (mean-field) theory, Eq. ([T]). Behavior in the 
different regions is as follows: a) stable, non-oscillatory 
(both eigenvalues real and negative), b) stable, oscillatory 
(eigenvalues complex, negative real part), c) unstable, os- 
cillatory (eigenvalues complex, positive real part), d) un- 
stable, non-oscillatory (both eigenvalues real and posi- 
tive). Note that there is a stable limit cycle in both c) 
and d). 

G {1, . . . , 4}, which are in turn functions of the pop- 
ulation vector n. The reactions of the Brusselator 
are given here alongside their corresponding transi- 
tion rate. 
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Note that different microscopic formulations of the 
Brusselator can be foimd in the literature (see e.g. 
[26, 27, 28]), with their differences pertaining only to 
the precise role of the non-dynamical substances. We 
here follow the formulation of [23|. 

The rates, T^ija), are derived according to the 
stoichiometric weights of each type of molecule on 
the left-hand side of each reaction. The first reac- 
tion effectively corresponds to the (spontaneous) cre- 
ation of Xi molecules. The rate of this reaction is 
proportional to the number of A molecules which 
we parametrize by the integer N and hence we use 
ri(n) ~ N . The spontaneous decay of Xi molecules, 
Xi —f 0, occurs with a rate proportional to ni and 
will be set to 12(11) — ni. There is no loss of general- 
ity by doing this since we can consider this reaction as 
setting the time scale of the model. With the choices 
made so far, the first two reactions in isolation ensure 
that the mean number of Xi molecules over time is 
TV (see for example The third reaction con- 

verts molecules of type Xi into type X2 and occurs 
at rate 13(11) = bni, where 6 is a parameter equal to 
the ratio of the number of B molecules to that of A 
molecules. The fourth reaction is of an auto-catalytic 
nature (with Xi being both a reactant and a product 
of the reaction) and will convert X2 into Xi at a rate 



4 



which is quadratic in ni and Hnear in n2. More pre- 
cisely, the rate of this fourth reaction is proportional 
to ni(ni — 1). However, we shall assume that ni » 1 
at all times and that T4(n) is well approximated by 
T4(n) = cN~'^n\n2- The factor of N~'^ ensures that 
c is of the same dimension as the parameter h. So 
long as h and c are independent of N , n2 also scales 
with N and hence N controls the total number of 
particles in the system. This justifies the use of the 
label system size for the parameter N. 

The stochastic time evolution of the system can be 
described by a master equation. To compactify no- 
tation we will encode the effect of an occurrence of 
reaction v G {1, ... ,4} on the system in a vector Vj^ 
describing the change of populations due to this reac- 
tion. The first component of v^, denotes the change 
in the number of Xi molecules due to the single oc- 
currence of reaction and the second component the 
change in the number of X2 molecules. For example, 
an occurrence of the first reaction increases ni by one 
while leaving n2 unchanged. In all we have, 

VI = (1,0), V2 = (-l,0), 

V3 = (-l,l), V4 = (l,-1). (5) 

The evolution of the time-dependent probability, 
^n(i), of finding the system in state n = (71,1,712) 
at time t is described by the master equation, 

, 4 

^Pn(i) = (^-(" - v,)P„_v„(0 - r,(n)P„(t)) , 

(6) 

subject to initial conditions Pn(^o)- The first term in 
the summation accounts for transitions of the system 
state from n— Vi, —^ n while the second term accounts 
for transitions away from the state n. 

While the problem of solving the master equation 
exactly is intractable, the mean-field behavior of the 
system, Eq. ([T]), may be recovered by multiplying 
the master equation on both sides by the population 
vector n followed by a summation over all possible 
configuration states. This leads to 

A(n)=^v.(T,(n)), (7) 

where the brackets (. . .) denote a time-dependent en- 
semble average over realizations of the stochastic dy- 
namics, i.e. (/(n)) = J2n -^nit) f (n) for any func- 
tion /(n) of the state-vector. We shall now also 
use the mean-field approximation {Ti,{n)) « T^((n)), 
which amounts to neglecting correlations by replac- 
ing {n\n2^ with (tt-i)^ ("2)- To simplify notation 
further we define the dimensionless variable x(t) = 
{n{t)) and dimensionless transition rates a,y(x) 
through the identity Na^{x.) = T^{Nx.) for all reac- 
tions. Finally, the rate of change of the mean-field 



concentration is 

X = A(x) = v,,ai,(x), (8) 

V 

where v runs over all reactions. In the Brusselator 
we have 

ai(x) = l, a2(x)=a;i, 
03 (x) = hx\, 04 (x) = cx\x2^ (9) 

and so we have the functions >li(x) = 1 — a;i(l + fe — 
CX1X2) and ^2(x) = x\(h — CX1X2) which agrees with 
the right-hand sides of the mass-action equations (see 
Eq. 111)). We shall sec that the definition of a^{x.) 
finds further use in the following section. 



III. STOCHASTIC EFFECTS IN THE 
FIXED-POINT PHASE 

One of the aims of this paper is to understand how 
the cycles generated by stochastic amplification, in 
the regime where a stable fixed point exists, behave 
as parameters change so that the fixed point becomes 
unstable and a limit cycle is born. Therefore in this 
section we examine the nature of these stochastic cy- 
cles by restricting the analysis to choices of the model 
parameters corresponding to points in phase b) of 

Fig. m 

In the previous section we showed that the average 
dynamics of a large ensemble of the finite Brusselator 
system, defined by the reaction dynamics of Eq. ([4|. 
will follow the deterministic path laid out by Eq. 
Individual realizations will, of course, follow random 
paths but simulations show these will move towards 
and subsequently stay close to the fixed point. In 
particular the random fluctuations away from a sta- 
ble fixed point solution are of relative order l/i/iV, 
as has been observed in similar analysis of other re- 
action systems [1, H, [13, HI] . The characteristics of 
the fluctuations can be studied analytically by means 
of an expansion of the master equation in the inverse 
system size ^] . This is a standard tool in the analy- 
sis of interacting particle systems, commonly referred 
to as van Kampen's system-size expansion, which has 
been applied to a number of different systems e.g. in 
[1, H, [TH, HI] , so we will not present the full details 
of the mathematical analysis here. Some of the inter- 
mediate steps are reported in the Appendix. The key 
idea is to write the particle populations ni and 712 of 
the finite system as 

m/N = xI+^i/Vn, 

n2/N = x*2+^2/Vn. (10) 

It follows that and £,2 are also random variables 
which represent the fluctuations of the dynamics of 
the finite system about the stationary solution of 
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the mean- field equations. We will frequently use the 
shorthand ^ = (i^i,^2)- A systematic expansion of 
the master equation ([6]) in powers of 7V~^/^ is then 
carried out following the lines of [1^. The leading 
order terms yield the deterministic mean-field equa- 
tions while the next-to-leading order terms give rise 
to a linear Fokker-Planck equation (see Eq. (jA-Sp ). 
describing the time evolution of the probability den- 
sity function of ^. The corresponding drift matrix is 
the Jacobian K* of the deterministic dynamics (as 
defined in Eq. The appearance of K* should 

not be surprising: Eq. (fTU]) has the same form used 
in the linear stability analysis of Section III Al with 
e = The diffusion matrix, D* , in the Fokker- 

Planck equation is given by Eq. (|A-10p evaluated at 
the fixed point. This yields, 

The Fokker-Planck equation (IA-8|) is equivalent to a 
Langevin system of the form 

m = K*m + m, (12) 

where i{t) = {fi{t), f2{t)) is bivariate Gaussian white 
noise of zero mean and with the following co- variance 
matrix indicating correlations between components: 

{Mt)fj(t')) - 2D*^5{t - t') e {1,2}. (13) 

Due to the linear character of Eq. (fT^ and given 
that the drift matrix is constant in time it is straight- 
forward to obtain analytical expressions for the power 

spectra Pi{u;) = and P2(^) = {Muj)\^). 

We have here written ^i(w) for the Fourier transform 
of the fluctuating variables ^i(t) {i — 1,2) with, 

/oo 
6(i)e-'-*dt. (14) 
-oo 

Following the steps of 0, H, [l^ one finds 

Pi(c^)= 2((l + 6)c^2^c2)l?-i(w), (15) 
P2{uj)= 2b{uj^ + l + b)V-\uj), (16) 

V{UJ)^ (c-Lj2)2 + (l+c- 6)2^2^ (17) 

As seen in Fig. [2]these spectra each show a maximum 
at a non-zero frequency, indicating amplified coher- 
ent oscillations due to the demographic noise. These 
analytical predictions compare well against simula- 
tions for different values of the model parameters b 
and c well inside the fixed-point phase. Numerical 
estimates for the power spectra are obtained through 
the repeated simulation of the microscopic chemical 
reactions using the Gillespie algorithm [26(. This is 
a widely used method to sample random paths from 
the solution to a master equation derived for Marko- 
vian particle systems. Only as the boundary of the 
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FIG. 2: (Color on-line) Power spectrum Pi{(^) = 

1^1 (w)!^^ of fluctuations in the concentrations of Xi 

molecules in the fixed-point phase (c = 1,6 = 
1.8, 1.85, 1.9, 1.95 from bottom to top at the maximum). 
Solid lines show results from the analytical theory, mark- 
ers are from stochastic simulations using Gillespie's algo- 
rithm (simulations are run up to t{ = 150, system size is 
A'^ — 10^, averages over 10* samples are taken). 

: » »b=i.9 : 
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FIG. 3: (Color on-line) Power spectrum Pi{u!) = 

|^i(a;)P^ of fluctuations in the concentrations of Xi 

molecules in the fixed-point phase near the onset of limit 
cycle behavior. Solid lines show results from the analyt- 
ical theory, markers are from stochastic simulations (run 
up to tf = 150, system size is A'' = 10^, averages over 10'* 
samples are taken). 

fixed point phase is approached (i.e. as 5 — > 2 from 
below for the fixed value of c = 1) do systematic de- 
viations between the theory and Gillespie simulations 
emerge visibly in Fig. [31 In particular, the power 
spectra from simulations begin to exhibit peaks at 
harmonics of the fundamental frequency given by the 
first maximum, which are not captured by the linear 
theory. These are the early precursors to the onset of 
limit cycles due to the stochastic broadening of the 
Hopf bifurcation. Further discussion of the effect of 
stochasticity on a Hopf bifurcation can be found in 

One of the main points of interest is to see what 
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happens as the boundary of the fixed point phase 
is reached. When b = 1.95, for example, the peak 
can be seen to reach a height of 3 x 10"^. In fact, 
it appears from Fig. [3] that as 6 ^ 2 the maxi- 
mum of Pi{uj) tends towards infinity. Stochastic ef- 
fects and resonant amplification of oscillations can 
hence become macroscopic (i.e. of the same order 
of magnitude as the mean-field dynamics) for sys- 
tems of very large system size, close to the transition 
into limit cycle behavior. This is an extreme case 
of stochastic amplification due to a resonance as can 
be seen from Eq. P7|) : if c is set equal to 1 and 
b = 2 — (5, the denominator vanishes at frequencies 
given by = [1 - ((5^/2)] ± i6. When b < 2, there 
is no zero for real lo, however as 6 — s- 2, the pole 
approaches the real axis and becomes real at w = 1 
when b — 2. The linear stability analysis can be ex- 
tended into region b > 2: for & = 2 + 6, there is an 
unstable spiral of period 2tt. For very small S, when 
the exponential growth can be neglected, this is a 
center of period 27r, which is the nascent limit cycle 
— not to be confused with a perturbation about the 
limit cycle to be discussed shortly. 

The next section will discuss general technical de- 
tails of how to characterize the stability of limit cy- 
cles in dynamical systems, and we will in particu- 
lar review elements of Floquet theory and Frenet co- 
moving frames. Both of these are standard tools used 
to study dynamical systems exhibiting limit cycles, 
and are as such not directly concerned with stochas- 
tic effects, but with perturbations and fluctuations 
about periodic attractors in general. We will return 
to stochastic systems in Section|Vl before conclusions 
will be drawn in Section IVTl 



IV. FLOQUET THEORY AND ROTATION 
INTO FRENET FRAME 

A. Limit cycles in the Brusselator system 

In the phases labeled by c) and d) in the phase dia- 
gram (Fig. [1]), the deterministic Brusselator system, 
as described by Eq. ([T]), exhibits limit-cycle behav- 
ior since we can eliminate both stable fixed points 
and unbounded trajectories. For positive initial con- 
ditions, Eq. ^ admits a stable periodic solution of 
period T. The period T will generally depend on 
the choice of model parameters b and c. It is found 
that, although T changes significantly with c, if we 
set c = 1 as we do in this paper, and change b only 
in a narrow band about 6=2, the period will change 
little from the value of 2tt found in the last section. 
Thus in what follows we will find that the angular 
frequency of the limit cycle remains close to w = 1. 

We will label limit cycle solutions by x(i) = 
{xi(t),X2(t)) in the following calculations and have 
x(t + T) = x.(t) for all times, t. In general the curve. 
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FIG. 4: (Color on-line) Illustration of the limit cycle solu- 
tion {xi{t),X2{t)) of the deterministic Brusselator system 
at fixed parameters 6 = 2.2, c = 1. 



x(t), cannot be calculated in closed form. However 
good estimates can be obtained via numerical inte- 
gration of Eq. 111). The geometrical shape of the 
limit cycle of the Brusselator model is illustrated for 
a fixed choice of the model parameters b and c in Fig. 

m 

In order to study stability, we now consider a dy- 
namical path beginning close to, but not on, the limit 
cycle, x.{t). If the limit cycle solution is stable then 
the difference between this path and the geometric 
curve of the limit cycle will decay as time progresses. 
Similarly to the expansion about a fixed point, we 
can write this difference as 

em = x(t) - x(t) (18) 

where, again, e expresses our anticipation that the 
deviation from the limit cycle is small. Expanding 
Eq. ([T]) in powers of e and letting e — > 0, one then 
finds that the time evolution of ^{t) takes on the lin- 
ear form, 

±m = Kitm, (19) 

where the matrix K{t) is found to take the specific 
form of Eq. (|A-9p for the Brusselator model. Study- 
ing the local stability of limit cycle solutions against 
perturbation is hence the analogue of studying the 
stability of fixed points as discussed above. The ele- 
ments of K{t) are given by Kij{t) — Kij(x{t)) which 
is simply the matrix, Kijijx.) — dap{ii) / dxj evaluated 
at the limit cycle. Therefore, due to the periodic na- 
ture of x(t), all elements of K{t) are periodic. 
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B. Floquet Theory 

An analytical tool to characterize the stability or 
otherwise of limit cycle solutions is Floquet theory — 
the mathematical theory of linear differential equa- 
tions with periodic coefficients (see |3C|, whose no- 
tation we will use). Since, in Eq. (jl9p . we have 
K{t + T) — K{t), Floquet theory is applicable. In 
our case, T is the period of the mean-field limit cycle 
under consideration. 

In essence Floquet theory states that, provided 
X{t) is a fundamental matrix of the system (fT9| . then 
there exists a non-singular constant matrix B such 
that 

X{t + T)=X{t)B (20) 
for all t. In addition one has 

detB = exp tvK{t)dt^ . (21) 

While the matrix B in general depends on the choice 
of the particular fundamental matrix X(t), its eigen- 
values (and determinant) do not. The eigenvalues of 
B are usually referred to as the Floquet multipliers 
of the system (fT9|). In the case of the Brusselator the 
matrices K{t), X{t) and B are 2x2 matrices and we 
denote the resulting Floquet multipliers by pi and p2 ■ 
Characteristic exponents pi and fi2 are then defined 
by setting pi = e^'^ for i £ {1, 2}. Further results of 
Floquet theory then concern the solutions of p^ . If 
yO is a characteristic multiplier for and p the cor- 
responding exponent then it can be shown that there 
exists a particular solution ^{t) of (|19p . such that 

i{t + T)=pm Vt. (22) 

One then finds that this solution can be expressed 
in terms of a periodic function p{t) (i.e. one with 
p{t + T) = p{t)) scaled by an exponential, 

m = e^*p(t). (23) 

General solutions of (fT5|) can therefore be written as 
a linear combination of functions of this form. For 
example in our two-dimensional system, 

|(t) = cie'^i*p(i)(i) + C2e'^^*p(2)(t), (24) 

with Ci , C2 constant coefficients determined by initial 
conditions. 

The Floquet analysis simplifies for the class of 
problems where the linear differential equations ([TO]) 
are derived from a dynamical system, x(t) = A(x), 
with a limit cycle x(t). In this case, it is easy to see 
by differentiation of the original equation of motion 
that the vector of velocities, x{t) = {xi{t),X2{t)), is 
a solution to (I19p . Since the velocity vector itself 
is a periodic function of time, we are therefore as- 
sured that one of the Floquet multipliers is equal to 



unity, pi = 1. That is, the corresponding exponent, 
pi, vanishes. This is a general result for all linear 
expansions about limit cycles arising from first-order 
equations. The remaining eigenvalue of B can then 
be determined using Eq. (PT|) and specifically for the 
Brusselator system we find that the corresponding 
Floquet exponent is given by 

1 /"^ 

P2 = 7t; / i-l-b + 2cxiit)x2{t)~cxi{t)^) dt. (25) 
^ Jo 

This integral can be evaluated numerically for any 
choice of the parameters 6, c which give rise to a limit 
cycle. For b = 2.2, c = 1 one finds p2 = -0.20225 
along with the already established observation that 
pi = 0. The corresponding functions p^^^(t) and 
p(^^ {t) are illustrated in Fig. [5] In fact for the Brus- 
selator the non-zero exponent is bound to be real and 
negative throughout phases c) and d) of the phase di- 
agram (Fig. [T]), i.e. throughout the limit cycle phase. 

In conclusion, we have established that one of the 
Floquet exponents of the system vanishes throughout 
this phase and that the remaining exponent assumes 
negative real values. The zero exponent is associated 
with perturbations in the longitudinal direction of the 
limit cycle; such perturbations are neither amplified 
nor reduced as the motion progresses. Perturbations 
in the transverse direction, by contrast, decay in time 
in the Brusselator system, rendering the limit cycle 
stable. Indeed, the multiplier p2 can be seen as char- 
acterizing a Poincare map of transverse motion. If the 
system is perturbed transversely by a small amount, 
(5, at time t = one may construct a Poincare map in 
the usual way [3l| : by forming the line perpendicular 
to the limit cycle which includes the point x.(t = 0). 
Then at every integer multiple nT of the period of 
the limit cycle the trajectory intersects the line at a 
distance P2<^ from the limit cycle. Since p2 < 1, this 
approaches the limit cycle with increasing n. 

C. Rotation into Frenet co-ordinates 

As seen above the Floquet exponents and periodic 
functions of the Brusselator system describe the re- 
laxation of perturbations in longitudinal and trans- 
verse directions. It is hence convenient to study the 
dynamics of deviations in co-ordinates defined along 
the tangential (longitudinal) and normal (transverse) 
directions to the limit cycle. This co-moviiig frame 
is generally referred to as the Frenet frame We 
will label the transverse coordinate by r and the lon- 
gitudinal one by s, as illustrated in Fig. [6l We 
follow the procedure of [13] and define (j){t) to be 
the angle between the x-axis and the line normal 
to the limit cycle at time t. Transformation of per- 
turbation displacement vectors from local Cartesian 
co-ordinates ^ — (^1,^2) into the co- moving Frenet 
frame q = (r, s) can then be thought of as a rotation 
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FIG. 5: (Color on-line) Periodic functions of the Floquet 
analysis of Eq. (I19p . The curves show each component 
of the vector- valued functions p'-^'(t) and p'^-'(t) respec- 
tively. The function associated with the vanishing expo- 
nent is displayed in the upper graph, while that associated 
with the negative real exponent is displayed in the lower 
graph. 




FIG. 6: Illustration of the rotation from Cartesian co- 
ordinates (^1,^2) to the Frenet frame, denoted by (r, s). 



by the angle 0(f), i.e. q(i) — J{t)^{t), where 

The angle (j){t) , in turn, can be expressed as a function 
of the velocity vector x.{t) at all points on the limit 
cycle trajectory. The transformation matrix is then 



w,^ _ J_ f X2it) -xi{t) 

''^'> v{t) \Mt) Mt) 



(27) 



where v{t) is the speed (magnitude of the velocity) 
given by v{t) = \/ (xi)^ + (a;2)^. The time derivative 



of deviations from the limit cycle expressed in the 
Frenet frame is then given as 

q(i) = J{t)i{t) + j{t)m 

= {Jit)K{t)J-\t)+R{t))ciit), (28) 

where R{t) — J{t)J^^{t) has been introduced and 
explicitly we have, 



R{t) 



X2{t)xi{t) ~ Xi{t)x2it) 



-1 

1 



(29) 



To simplify notation wc will use K'{t) = 
J{t)K{t)J-\t) and K^°\t) = K'{t) + R{t). As seen 
in [22 the upper-right element of if vanishes iden- 
tically, i.e. we have Klg{t) + Rrs{t) = for all times. 
Hence, the motion of the first co-ordinate r{t) in the 
Frenet frame decouples from the second. More pre- 
cisely, Eq. (fT^ takes the form 



d_ 

dt 



r{t) 
s{t) 



Klf{t) 
Klf(t) Klf{t) 



r[t) 
s{t) 



(30) 



after rotation into co-moving co-ordinates. The non- 
trivial elements of i<f *°* can be computed explicitly as 
functions of the limit cycle trajectory and are given 

by 



1 



2 {2^1-^22 -^2 



xIKii 



- a::ia;2(-ftri2 + ^^21)}, 
+ 2icii:2{Kii- K22)}, 



KIT = v/v. 



(31) 

(32) 
(33) 



D. Re-scaled Frenet frame 

Further simplification can be achieved by re-scaling 
the co-ordinates of the Frenet frame, after rotation, 
by the velocity of the limit cycle. More precisely we 
make the transformation x{t) = 'i{t)/v{t), so that 

v{t) 

We will denote the individual components by x(t) = 
{p{t),a{t)). One finds that the perturbative displace- 
ment is described by 



xit) 



^/)x(t)^i*°'xW, (35) 



where / denotes the identity matrix. Combining Eqs. 
((501) and Eq. ^ takes the simple form. 



A { Pit) 

dt V <t) 



LIT it) M (pit) 



(36) 
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variables n i-^ ^{t) defined by 



-0.5 



FIG. 7: (Color on-line) Periodic functions p^p\t) and 
pl?\t) associated with the non-trivial Floquet multiplier 
p2 of Eq. (|36|l . While a perturbation p decays in a non- 
linear fashion, the phase, cr, oscillates. 



In particular the lower right element of L*°* vanishes 
identically. The remaining elements are given by 



-2^,MKi2 + K2i)}, (37) 

C;(t) = l{(5^_5;i)(Xi2 + if2i) 

+ 25;i5;2(i^ii-if22)}. (38) 

These relations are valid for general systems of first- 
order ordinary differential equations with two degrees 
of freedom which exhibit limit cycle solutions. The 
vanishing elements in Eq. (|36|) guarantee the exis- 
tence of the constant solution x(t) = (0, ctq) (with ctq 
a real- valued constant). Hence we find that a pertur- 
bation in (T is not only periodic, as indicated by the 
trivial Floquet exponent, but it is in fact constant. 
Transverse perturbations relax in a non-trivial way 
and for completeness we show the associated periodic 
functions, Pp'\t) and pi^'(t) in Figure [7l The super- 
script here indicates that we refer to the non-trivial 
Floquet multiplier p2- It should be noted that there 
is a secondary oscillatory effect on cr as p decays back 
to the limit cycle. Since this velocity-scaled rotation 
yields a simpler linear theory than rotation alone, we 
will only use the co-ordinates (p, cr) in the discussions 
that follow. 



V. STOCHASTIC EFFECTS IN THE LIMIT 
CYCLE PHASE 

A. System size expansion and analytical 
predictions 

We will now return to the stochastic system defined 
by the chemical reactions (HJ) or, equivalently, by the 
master equation In this section we will apply van 
Kampen's system-size expansion to the case where we 
have a limit cycle, (a;i(i), X2(i)), in the mean-field. 
The starting point is the transformation of random 



n = iVx(t) + \/Ni{t) 



(39) 



The only difference compared to the transformation 
(jlOp we applied in the fixed-point phase is the time 
dependence of the first term on the right-hand side. 
Apart from this minor complication the algebraic 
steps necessary to carry out the expansion of the mas- 
ter equation are mostly unchanged, see the Appendix 
for details. As before, the evolution of the stochastic 
fluctuations is described by a linear Langevin equa- 
tion. 



m = K{t)£,+m, 



(40) 



where f (i) — (/i(i), f2{t)) is bivariate Gaussian white 
noise with zero mean with correlations given by, 
{Mt)fj{t')) = 2D^j{t)S{t - t') [ll, 22j. The forms of 
the matrices Kit) and D{t) are given in the Appendix 
(Eqs. (|A-9|) and (IA-10|) . respectively). Because of the 
periodicity of x(f;), all elements of Kit) and D{t) are 
periodic with the period of the limit cycle, T. 

The linearity of these Langevin equations allows 
us to make further analytical progress. To this end, 
it is convenient to study the stochastic dynamics in 
the velocity-scaled Frenet frame as introduced above. 
Upon performing a rotation into the co-ordinate sys- 
tem spanned by (p, cr), Eq. (jlU]) takes the form 



d_ / pit) 
At \ a{t) 



'(t) 



T tot 



Pit) 
a{t) 



(41) 



with {C^{t)) =0, i = 1,2 and 



mW)) = 2[K{t)D{t)h'{t)],,5{t - t'). (42) 

Hence the correlations of the noise components in the 
co-moving frame are described by the matrix H{t) = 
Ki1:,)D{t)K^{t), where K{t) is defined by Eq. ([51]) . 

Since the mean values of both p and cr vanish, 
one can write the variances of these variables as 
Vpp{t) = {p^{t)) and Vaa{t) = (cr^(i))- First-order 
ordinary differential equations can be derived [2^ for 
these quantities, and they take the form 



= 2L\°;Vpp{t) + 2Hpp{t), 

'tot 



(43) 



L':'^Vpp{t) + L*;;Kp(i) + 2H,pit)m 

(45) 



2L'°'pV,p{t) + 2H„„{t). 



Vppit) 

Kit) 

Kit) 



Solving Eqs ([43 |1 - ([45|) sequentially, a closed form for 
Va-crit) may be found and evaluated numerically. Re- 
sults are shown in Fig. [S] Note that the variance of 
cr, evaluated at integer multiples of the time period, 
T, increases linearly. This increase without bound 
means that the linear approximation within which we 
derived our theoretical results can be expected to be 
valid only at sufficiently short times. More precisely 



'0 




FIG. 8: (Color on-line) Variance Vcrcr{t) of longitudinal 
fluctuations as obtained from solving Eqs. (|43p -(|45 p . 
Model parameters are fixed to b = 2.2, c = 1. The solid 
line is calculated from the theory while the markers are 
obtained from simulations. System size is A'' = 10^ , aver- 
ages over 10000 runs are taken. 



the first-order van Kampen expansion is accurate pro- 
vided (j{t)/\/N is small compared to the components 
of the hmit cycle solution, x(t). The time scale on 
which longitudinal fluctuations remain small enough 
for the linear theory to be valid will hence increase 
as the system size is increased, and will diverge as 
N ^ oo. 

The Langevin equation for radial fluctuations com- 
pletely decouples from that of a. From Eq. ((4T|) we 
have. 



\t)p(t)+(^,{t), 



(46) 



which is readily integrated, given an initial condition 
p{to) = Po- 



pit) - pQ<^(t,to) 



^t,t')C,{t')dt', (47) 



where we have used the definition 

Ht,t') = e^p(^l\xit")dt"y (48) 

We can now evaluate the average temporal correla- 
tions of radial fluctuations for r > 0: 

(p{t + T)p{t)) - po'^it + r, i)$2(t, to) 



Jto 

= 2^t + T,t) [ <i>\t,t')Hpp{t')dt', (49) 

Jto 

where we have used the identity (t){t,t')(t){t' ,t") = 
(j){t,t"), valid for all t,t',t". We now set the initial- 
ization time to the inflnite past [to — > — oo) so that 
the initial condition itself is forgotten. After invok- 
ing the periodicity of Lpp(t) and Hpp{t) and making 
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FIG. 9: (Color on-line) Power spectra Pp{i^) and Pr{io) of 
transverse fluctuations in the velocity scaled Frenet frame 
and ordinary Frenet frame for the Brusselator at c = 1 
and b = 2.2. The solid pair of curves shows results for 
Pp{u}) while the dashed- line pair shows those without the 
rescaling of the Frenet frame. The monotonic curves in 
each case correspond to Lorentzians obtained by replacing 
V-°p{t),Hpp{t) and Kl°^,Grr{t) by their time averages in 
Eq. 



a suitable change of the integration variable one finds 



{pit + T)p(t)) 



2$(t -I- r, t) 



T 

XI <i>^it + T,t + t')Hppit + t')dt', {50) 



where the pre-factor is due to an infinite summation 
of powers of e^^p , resulting from the observation that 
^{t + T,t) = e'^f, where pp is the non- vanishing Flo- 
quet exponent of the system. Eq. ([SU]) confirms that 
the two-time correlation {p{t + T)p{t)) is periodic in 
t. Averaging over t yiel ds the (time-averaged) auto- 
correlation function [l7l [s^ , 

C{t)^^ j\p{t)p[t + T))dt. (51) 

The power spectrum of p{t) is now obtained as the in- 
verse Fourier transform of this function C(t). A sim- 
ilar procedure can be implemented to compute the 
power spectrum of r(t), the only necessary replace- 
ments being Hpp{t) Grr{t) = [J{t)D{t)j'^ {t)]rr 

andL^°;(t)-.K^°*(i). 

The resulting spectra of transverse fluctuations are 
plotted in Fig. El both in the rescaled Frenet frame 
(the power spectrum is then denoted by Pp{uj)), and 
in the ordinary Frenet frame, see the curve referred 
to as Pr{oj)- Peaks of a finite width are found at 
multiples of the frequency of the limit cycle and it 
can be seen that the structure of these peaks is de- 
pendent only on the parameters b and c and not on 
system size, N. For comparison we also show the 



FIG. 10: Illustration of how data is obtained from Gille- 
spie simulations. At any time t the algorithm generates 
a point {n\{t) / N ,n2{t) / N) in the a;i, X2-plane. The vec- 
tor components of the vector ^(t) are then obtained as 
= N^'^ [ni{t)/N -xi{t)] and similarly for £,2{t). As 
before, x(t) is the point marked by the limit cycle tra- 
jectory at time t. Subsequently, {(t) is converted into 
co-ordinates = (p(^)i'^(^)) by the transformation 

X{t) = A(t)^(f), with ^{t) as defined in Eq. 



power spectra one would obtain by replacing the pe- 
riodic matrices L^°\t), K^°\t), H{t) and G{t) by 
their mean values (e.g. H = T^^ H{t)dt) in 
Eq. ([50]) . In this case the power spectra reduce to 
Lorentzian curves which we plot in Fig. [9] along with 
the full results. It would seem a plausible conclu- 
sion that the time-dependence of the drift and dif- 
fusion matrices (_L'°*(i) and H{t), respectively) con- 
tributes additively to the Lorentzian shape that we 
would expect from a Langevin equation with time- 
independent drift and diffusion matrices. 



B. Test against numerical simulations 

The above theoretical prediction for the power 
spectrum of transverse fluctuations about the mean 
field behavior can be verified in simulations of finite 
systems using the well known Gillespie algorithm [2^ 
to simulate the reaction system of the Brusselator 
model. This method generates continuous-time real- 
izations {ni{t),n2{t)) of the multi-particle stochastic 
process described by the master equation Note 
that ni{t) and n2{t) are integer valued at all times. 
The overall system size, N, is a control parameter in 
these simulations, as are the reactivity parameters b 
and c. 

To make contact with the above theoretical analy- 
sis we begin by considering the transformation n i-^ ^ 
given in Eq. i.e. we have, 



n2{t)/N = X2{t)+Ut)/N^'^- 



(52) 



The quantities p and a are then obtained by 
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FIG. 11: (Color on-line) Power spectrum Ppiuj) of trans- 
verse fluctuations about the limit cycle trajectory. Model 
parameters are c = 1,& = 2.2. The solid curve results 
from theoretical analysis while filled symbols are from 
Gillespie simulations run up to tt = 25 and connected 
open markers are from simulations run up to tf — 200. 
Simulations interpreted using the rotation method as ex- 
plained in Section FV Bl Averages over 10000 samples are 
taken, system size is TV = 10^. 



performing the transformation {p{t),a{t))'^ = 
A(i)(6(i),?2(t))'^, with A(t) as defined in Eq. 
This is illustrated in Fig. [TUl 

The linear theory we have developed can be ex- 
pected to be accurate at most in a regime where the 
second terms on the right-hand sides of Eq. ([52]) are 
both small corrections to the first terms. Equivalently 
we require that p{t) and a{t) both remain small com- 
pared to N^^'^. However, as we have already discussed 
and demonstrated in Fig. [51 the variance of a{t) 
grows linearly in time (modulo periodic variations). 
Hence the results based on the van Kampen expan- 
sion are expected to be accurate at most on time 
scales tf — 0{N'^). Performing simulations at system 
size N and run up to time scales tf ^ N one observes 
good agreement with the predictions of the linear the- 
ory as illustrated by Fig. Ill|, where we show data for 
systems of size N = 10^ run up to tf = 25. How- 
ever, when we consider simulations run up to larger 
times (e.g. tf = 200), the power spectrum shows sys- 
tematic deviations from the theoretical curve. This 
effect can be accounted for by the reasonably steep 
average slope found in Fig. [Sj Extrapolating Fig. 
[HI to larger times one expects {a{t = 200)^) w 4000, 

(cr(t = 200)/ V^)2\ w 0.2 for iV = 10^ so 

that the second terms on the right-hand sides of Eq. 
((52|) can no longer be thought of as small compared to 
the first. The next section will discuss an alternative 
set of measurements that may be taken from Gille- 
spie simulations. These will be shown to successfully 
tackle this defect and yield a good match with the 
prediction illustrated in Fig. [5] also on longer time 
scales. 
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FIG. 12: The modified procedure of identifying deviations 
from the limit cycle trajectory. The point of reference 
is now the point on the limit cycle trajectory, which is 
geometrically closest to the point in the Cartesian plane 
obtained in Gillespie simulations. 



C. Modified method of comparison 

The approach that we took in the previous section 
to interpreting time series from the GiUespie simu- 
lations compares extremely well with the theory we 
have developed in terms of the variance of a. How- 
ever, except for simulations performed on very short 
time scales, we find that this is not the case for the 
power spectrum of p. We will therefore pursue an al- 
ternative way of comparing data from Gillespie sim- 
ulations with the power spectra obtained from the 
theory. Specifically, the proposal put forward in this 
section is to use a transformation n i-^ k of the form 

n = N^{S) + %/iVK(n, S), (53) 

where 5' is a random variable chosen as 

S = arg min^, |7Vx(t') - n|. (54) 

Here arg min|Afx(t') — n| denotes the value t' min- 
imizing the quantity \N^[t') — n|, i.e. S is chosen 
such that x(S') is the point on the limit cycle trajec- 
tory with the minimum distance to the point n/7V 
obtained from the Gillespie simulation (| • | refers to 
the Euclidean norm). By construction the vector k is 
then perpendicular to the velocity of the limit cycle: 
/«(n, S').x(5') = 0. In this way we are able to define a 
Frenet frame directly in terms of the variable n which 
appears in the original master equation, rather than 
constructing it as a rotation from continuous Carte- 
sian coordinates. 

Although our intention is to use this formulation 
to reinterpret the simulation data, we first have to 
show how the above construction recovers the results 
obtained through the van Kampen system-size expan- 
sion. Since S* is a function of the stochastic variable 
n, we may define its mean value as 

(5)^ = ^5(n)P„(i), (55) 

n 



where Pn{t) is the solution of the master equation. 
Fluctuations about this mean value are defined by 

S{n) = {S)t + ^a{n,t). (56) 

Neglecting terms of order 1/N, as before, we may 
Taylor expand about the point {S)t on the curve to 
obtain 

nS)=n{S)t) + ^a{n,t)±{t). (57) 

Using Eqs. ([53]) and ([57]) we have that 
n = Nx{{S)t) + y/N {cr(n, t)±{t) + /«(n, S)} . (58) 

To the order we are working, we may replace n and 
S in the curly bracket by Nx.{{S)t) and {S)t respec- 
tively, by using Eqs. ([5^)1 and (|58p . Furthermore 
from Eq. fSS]) we have n — Nx{{S)t) to leading or- 
der. Since the limit cycle is defined by x, the correct 
mean-field equations are only recovered if the identi- 
fication {S)t = t is made. These considerations lead 
to Eq. ((55)) being written as 

n ^ NTi{t) + y/N {a{t)±(t) + K{t)} . (59) 

Comparing Eq. (I59p with the usual starting point for 
the van Kampen expansion, Eq. pop . we see that 
for these two approaches to agree it must be the case 
that 

^{t) = a{t)±{t) + K{t). (60) 

This is so, since to leading order the condition 
/t(n, S').x(S') — becomes K{t).x{t) = 0. Therefore 
from Eq. ([5U)) we require that a{t) = $,{t) .x.(t) / v'^ (t) . 
But from q{t) = J{t)^{t) and Eq. ^ it follows 
that s{t) = ^{t).ii{t)/v{t) and from Eq. ^ that 
a{t) = s{t)/v{t). This identifies <j as the scaled, lon- 
gitudinal component introduced in Section HV Dl By 
construction k = ^ — ctx must be the transverse com- 
ponent. In this way, we recover the van Kampen 
ansatz used previously. 

The Gillespie algorithm generates a time-series for 
n which we can use together with knowledge of the 
equation for the limit cycle to determine S, from 
Eq. ([54]) . This can then be used (i) to determine 
tr, which by Eq. ([55]) and the identification {S)t — t 
is y/N{S — t), (ii) to determine k from Eq. ([55)) . 
and (iii) to obtain re-scaled transverse fluctuations 
p{t) = 'ip{t)\K{t)\/v{S). The function ilj{t) is present 
to ensure the correct sign; we choose ip{t) = —1 when- 
ever n is located inside the limit cycle trajectory, and 
V'(i) = 1 if n hes outside the area surrounded by the 
limit cycle. It should be noted that now time is in- 
troduced only from the Gillespie algorithm, and not 
through Eq. (j53p . Using this methodology, we may 
determine the power spectra from repeated Gillespie 
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FIG. 13: (Color on-line) Power spectrum Pp{u->) of trans- 
verse fluctuations about the limit cycle trajectory. Model 
parameters are c = 1,6 = 2.2. The solid curve results 
from theoretical analysis while filled symbols are from 
Gillespie simulations run up to t{ — 25 and connected 
open markers are from simulations run up to tf = 200. 
Simulations are interpreted using the projection method 
as explained in Section fV CI Averages over 10000 samples 
are taken, system size is A'^ = 10^. 



simulations. The results for the variance of a are 
indistinguishable from those shown in Figure [S] and 
so we do not show them. However, in Fig. [T3] we 
compare the power spectrum estimated by measur- 
ing p{t) = ^p{t)\K{t)\/v{S) from simulations with the 
analytical results obtained in Section IV Al (see Eq. 

As seen in the figure, we find good quantitative 
agreement also on time scales on which the method 
discussed in the previous section failed to reproduce 
the theoretical curve. 

The findings of the section may be summarized by 
postulating three different temporal regimes. The 
first regime is defined for times which are sufficiently 
short that a/VN can be thought of as small com- 
pared with the size of the limit cycle. Thus Eq. ([5^ 
may be used, and the modified approach based on Eq. 
(f5S|) need not be used. The results for t{ = 25 exem- 
plify this regime. In the second regime, a j^pN is now 
sufficiently large that the use of the naive expression 
([5^ leads to a disagreement between simulations and 
the theoretical curve. This has been discussed above 
and is exemplified by the results for if = 200. Fi- 
nally, at longer times a will start to probe the peri- 
odic structure of the limit cycle, and a different type 
of behavior will occur. We have not explored this lat- 
ter regime in the present paper, but we will discuss 
it again in Section EH 

We end with a further illustration of the differ- 
ence between the two methods of extracting trans- 
verse fluctuations from simulation data. We show an 
example of a single run in Fig. [H] and depict the 
resulting two time series /9(i) produced using the two 
different methods. Initially the two time series agree 
well with one another. On longer time scales, how- 




FIG. 14: (Color on-line) Single time series of transverse 
fluctuations, p(t), illustrating the difference between the 
two measurement methods. The black curve shows data 
obtained using the rotation method (see Section fV Bp . the 
red curve shows result from the projection method (see 
Section FV C\ . Data is obtained from a single Gillespie run 
at Af = 10^ (6 = 2.2, c = 1). 



ever, systematic deviations are observed. We here 
note that the shape of time series as shown in Fig. 
[T3]can vary considerably between different runs of the 
stochastic Gillespie simulations. For reasons of clar- 
ity the realization shown in the figure is one where the 
deviation between the time series generated using the 
two methods is reasonably pronounced at large times. 
In other runs the discrepancy was smaller. 



VI. DISCUSSION 

In summary we have carried out an analysis of the 
effects of internal fluctuations found in finite systems 
with a large system size, N . We have focused on 
two-dimensional systems and we have used the Brus- 
selator model as a toy example which we have dis- 
cussed in two regimes. In the case where the mean- 
field dynamics approaches a stable fixed-point behav- 
ior with complex eigenvalues, we found the expected 
sustained oscillations driven by the stochasticity of 
the discrete particle dynamics. The power spectra of 
these oscillations can be obtained analytically via an 
expansion in the inverse system-size about the time- 
independent solution. This is similar to work carried 
out in the context of other models with fixed-point 
behavior 0,0,&[3. 

One of the aims of the present work was to extend 
these analytical tools to the case in which the dynam- 
ics exhibit a periodic solution on the mean-field level 
for a range of parameter values. The above van Kam- 
pen expansion in the inverse system size can be car- 
ried out as before. However, now that we are expand- 
ing about a curve, rather than a point, we must make 
a choice for the point on the limit cycle about which 
we expand. This scheme may be left undetermined 
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for carrying out the system-size expansion. However 
when carrying out simulations, more care has to be 
exercised when picking the point that one expands 
about. At short times, the time since the start of the 
simulation can be used to determine the point along 
the limit cycle. However, at longer times this will 
typically not coincide with the nearest point to the 
quantity n/N, due to the diffusion in the direction 
tangential to the limit cycle. We have given a pre- 
scription for carrying out simulations in the second 
temporal regime, and shown how we recover agree- 
ment between simulations and the van Kampen ex- 
pansion using this method. As we pointed out, there 
is a third regime where a/^/N becomes of the order of 
the period T, when further modifications will have to 
be introduced. Eventually, on very long timescales, 
the diffusive behavior can be studied by using N~'^ 
rather than N^^/"^ as an expansion parameter (25j |. 
It would be interesting to extend the work we have 
presented here to these longer times. 

Another objective was to understand the relation 
between the cycles due to the stochastic amplification 
in the fixed point phase and the limit cycles in the 
phase where the fixed point becomes unstable. We 
have shown how, for the example of the Brusselator, 
the former become the latter as one passes through 
the phase boundary. In the fixed point phase the fluc- 
tuations are amplified by a resonance which may be 
described by a pole in the complex frequency plane. 
As the phase boundary is approached the pole mi- 
grates towards the real axis, reaching it when the 
boundary is crossed and so turning the resonance 
into a limit cycle. Although this has been illustrated 
in the case of the Brusselator, we expect this phe- 
nomenon to be generic, and that it may be applied 
to the various systems mentioned in the Introduction. 

The output of the Gillespie algorithm is a time se- 
ries similar to that which is found in data obtained 
from real systems. We therefore expect that the 
methods we have applied in this paper will be ap- 
plicable to real data. We hope that this will lead to 
further insights when applied to the many other sys- 
tems which have a variety of stable attractors and 
which are subject to intrinsic noise. 
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APPENDIX: VAN KAMPEN SYSTEM SIZE 
EXPANSION 

In this appendix we briefly sketch some of the 
mathematical steps involved in carrying out the 



system-size expansion for two-dimensional chemical 
systems, of which the Brusselator is an example. The 
starting point is the master equation ([6]), 



At 



(T,(n- v,)P„_v„(i) - T,(n)P„(t)) , 

(A-1) 



which, subject to an initial condition Pnito), governs 
the temporal evolution of the probability distribution 
describing the statistics of the microscopic dynam- 
ics as defined by the reactions ([1]). We determine in 
Section |V] that we may use the following mapping 
between the population vector n and a continuous 
fluctuation, about the mean- field concentration 
as follows, 



(A-2) 



n = iVx(t) + V A^KO 



Following the procedure of |25| we now formulate the 
problem in terms the probability distribution H(^, i), 
describing the statistics of the stochastic process ^{t). 
Since ^(t) is a linear transformation of n we have 
that n(^,t) cx Pn{t). Hence we may directly sub- 
stitute n(^,t) into the master equation (|A-ip . We 
firstly note that the derivative with respect to time 
in (|A-ip is taken at constant n so that we have 
di/dt = -7Vi/2dx/dt. This leads to 



At 



aH(i,i) dxi{t) 



-N 



dt d^i 
./^ dUitt) dx2{t) 



dt 



d6 



dt 



(A-3) 



Next, we write the right-hand side of Eq. (|A-1|) in 
terms of H(^, t) and find that 



dt 



•^n^^'*) .5:[cxp(-A.-v....v,) 

X \jVa^{x{t) + N-^/^K)U{^,t) 



- 1 



(A-4) 



where we have also used the re-scaled reaction rates 
Cj/ as introduced earlier, i.e. 

T^{N^{t) + = Na^ (x(t) + . (A-5) 

We have exploited the continuous nature of ^ by using 
the shift operator: the exponential in the differential 

operator, = used in Eq. (jA-4|) . has 

the effect of shifting the argument of the subsequent 
functions by the vector —N^^^'^w^. Explicitly, for 
any smooth function, -F(^), 

exp(-iV-i/2v,-V^)F(0 = F{i-N-^/^^,). (A-6) 

Both the exponential and the rate functions, aj/(x), in 
Eq. (jA-4[) may be expanded as polynomials in N~^^^. 
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We then take the formal hmit of ^ oo keeping only 
terms of the two highest orders in N. One may then 
equate coefhcients of powers of N between the left 
and right-hand sides of the master equation (given 
by l|A-3|) and (|A-4p respectively). 

To leading order one consistently recovers the 
mean- field equations ([7]), i.e. 

^x(t)=A(x(^)), (A-7) 

with A(x) — ^j^v^a,/(x). Expanding to next order 
one finds 

+ EA.w|-|-n(tO, (A-8) 

which is a linear Fokker-Planck equation for the dis- 
tribution n(^,f). The elements of the drift ma- 
trix K(t) are given by Kij{t) — -^Ai{-x.{t)) and 



those in the diffusion matrix D{t) are Dij{t) = 

For the example of the Brusselator that we use, 
these matrices are explicitly given by, 

and 

^/ N _1/ l + a;i(l + 6 + CX1X2) -xi{b + 0x1X2) \ 
2\ —Xi{h CX1X2) Xi{b + cxiX2) J 

(A-10) 

Eq. (jA-8[) may be considered directly as a stochas- 
tic process ^(t) as described by a Langevin equation. 
This is discussed in the main text for two cases. In 
the case of a globally stable fixed point, x = x*, if(x) 
and D{x.) become K* and D* . In the case of a limit 
cycle, X = x(t), the drift and diffusion matrices are 
functions of time, K{t) and D{t), which naturally 
pick up the periodicity of this solution. 
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